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FINITE DIFFERENCE METHODS FOR THE INFINITY 
LAPLACE AND p-LAPLACE EQUATIONS 

ADAM M. OBERMAN 



Abstract. We build convergent discretizations and semi-implicit solvers for 
the Infinity Laplacian and the game theoretical p-Laplacian. The discretiza- 
^ tions simplify and generalize earlier ones. We prove convergence of the solution 

^ ^ of the Wide Stencil finite difference schemes to the unique viscosity solution 

^^ of the underlying equation. We build a semi-implicit solver, which solves the 

^^ Laplace equation as each step. It is fast in the sense that the number of itera- 

^^ tions is independent of the problem size. This is an improvement over previous 

^^ explicit solvers, which are slow due to the CFL-condition. 

< 

2 1. Introduction 

• 

f^ The Infinity Laplacian equation is at the interface of the fields of analysis, non- 

C^ linear elliptic Partial Differential Equations (PDEs), and probabilistic games. It 

G was first studied in the late 1960s by the Swedish mathematician Gunnar Aronsson 

I— I [Aro67, Aro68, Aro84], motivated by classical analysis problem of building Lips- 

chitz extensions of a given function. Aronsson found non-classical solutions, but a 

rigorous theory of weak solutions was not yet available. It took a few decades until 

QQ analytical tools were developed to study the equation rigorously, and computational 

t — ■ tools were developed which made numerical solution of the equation possible. 

CN In the last decade, PDE theorists established existence and uniqueness, and 

^'l regularity results. The theory of viscosity solutions [CIL92] is the appropriate 

t"^ one for studying weak solutions to the PDE. But the general uniqueness theory 

^^ did not apply to this very degenerate equation, so proving uniqueness required a 

^^ new approach. The first uniqueness result was due to Jensen [Jen93], followed by 

• • a different proof by Barles and Busca [BBOl]. Later, the connection with finite 

^^ difference equation was exploited by Armstrong and Smart, and they were to give 

S> a short uniqueness proof for the PDE [AS 10] . The differentiability of solutions 

%^ remained an open question for some time. The first result was obtained by [Sav05] in 

^ two dimensions, followed by [ES08], and [ESI lb] and [ESI la] in general dimensions. 

The first convergent difference scheme was presented one of the authors [Obe05] . 
A numerical scheme using the extension property can be found in LeGruyer [LG07]. 
Two different numerical methods were derived in [ESI la], one adapted the mono- 
tone scheme in [Obe05] to the standard Infinity Laplacian (which is homogeneous 
degree two in Vii), the second, quite different, used the variational structure of a 
regularized PDE. 
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Earlier work by LeGruyer [LGA98] proved uniqueness for a related finite dif- 
ference equation. The proof is a generalization of the uniqueness proof for linear 
elliptic finite difference schemes from [MW53]. 

A group of probabilists, Peres-Shramm-Sheffield- Wilson [PSSW09] studying a 
randomized version of a marble game called Hex found a connection with the Infinity 
Laplacian equation. This connection gives an interpretation of the equation as a 
two player random game. This is related to work by Kohn-Serfaty [KS06] who found 
an interpretation of the equation for motion of level sets by mean curvature, [OS88] 
[ES91] as a deterministic two player game. This equation was also interpreted as a 
nonlinear average [Obe04] for the purpose of finite difference schemes. Deterministic 
game interpretations for more general PDEs followed in [KSIO]. The connection 
between these various game interpretations was further studied in [EvaOT]. 

The rich connection between games, finite difference schemes, and nonlinear 
elliptic PDEs is now much better understood. There have been a number of works 
in this area, in particular on the game theoretical p-Laplacian. The probabilistic 
games interpretation can be found in [PS08] (see also [MPRIO]). Related works 
include biased games which corresponds to a gradient term [ASS 11] 

This article will further exploit the connection between games, finite differ- 
ence schemes, and nonlinear elliptic PDE, by building convergent finite difference 
schemes which are consistent with the game interpretation. The existence and 
uniqueness results are now established, efficient numerical solution of the equation 
remains a challenge. The original convergent scheme proposed in [Obe05] con- 
verged, but is not efficient: as the grid size grows, so does the number of iterations 
required to find the solution. This article improves and simplifies the original dis- 
cretization, and also finds fast solution methods. It also generalizes the scheme and 
the solvers to the game-theoretical p-Laplacian, which is a convex combination of 
the Laplacian and the Infinity Laplacian. 



1.1. Introduction to numerical methods for degenerate elliptic PDEs. 

There are two major challenges in building numerical solvers for nonlinear and 
degenerate elliptic Partial Differential Equations (PDEs). The first challenge is to 
build convergent approximations, usually by finite difference schemes. The second 
challenge is to build efficient solvers. 

The approximation theory developed by Barles and Souganidis [BS91] provides 
criteria for the convergence of approximation schemes: monotone, consistent, and 
stable schemes converge to the unique viscosity solution of a degenerate elliptic 
equation. But this work does not indicate how to build such schemes, or how 
to produce fast solvers for the schemes. It is not obvious how to ensure that 
schemes satisfy the comparison principle. The class of schemes which for which 
this property holds was identified in [Obe06], and were called elliptic^ by analogy 
with the structure condition for the PDE. 

An important distinction for this class of equations is between first order (Hamilton- 
Jacobi) equations, and the second order (nonlinear elliptic) case. The theory of 
viscosity solutions [CIL92] covers both cases, but the numerical methods are quite 
different. In the first order case, where the method of characteristics is available, 
there are some exact solutions formulas (e.g. Hopf-Lax) and there is a connection 
with one dimensional conservation laws [Eva98] . The second order case has more in 
common with divergence-structure elliptic equations, but because of the degeneracy 
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or nonlinearity, many of the tools from the divergence-structure case (e.g. finite 
elements, multi grid solvers) have not been successfully applied. 

In the first order case, there is much more work on discretizations and fast solvers. 
For Hamilton-Jacobi equations, which are first order nonlinear PDEs, monotonicity 
is necessary for convergence. Early numerical papers studied explicit schemes for 
time-dependent equations on uniform grids [CL84, Sou85]. These schemes have 
been extended to higher accuracy schemes, which include second order convergent 
methods, the central schemes [LTOO], as well as higher order interpolation methods, 
the ENO schemes [0S91]. Semi-Langrangian schemes take advantage of the method 
of characteristics to prove convergence [FF02]. These have been extended to the 
case of differential games [BFS99]. Two classes of fast solvers have been developed, 
fast marching [Set99], and fast sweeping [TCOZ03], The fast marching and fast 
sweeping methods give fast solution method for first order equations: both take 
advantage of the method of characteristics, which is not available in the second 
order case. 

There is much less work in the second order degenerate elliptic equations. The 
equation for motion by mean curvature [OS88], [ES91] has been extensively stud- 
ied. There is an enormous literature on this equation, but we just closely related 
references. The connection with games was already mentioned above. Numerical 
schemes include [CFF06] and [CFFIO] . In the case of motion by mean curvature, 
the equation is time-dependent, so a fast solver would allow larger time steps. For 
this equation, a semi-implicit solver has been built by Smereka [Sme03]. The idea 
from the Smereka paper will be adapted in this work to build fast solvers for Infinity 
Laplace. Another equation in this class is the Hamilton- Jacobi-Bellman equations, 
for the value function of a stochastic control problem. Applications include port- 
folio optimization and option pricing in mathematical finance. Numerical works 
include the early paper [LM80] and [CF95], and a paper on fast solvers [BOZ04]. 

For uniformly elliptic PDEs, monotone schemes are not necessary for convergence 
(for example most higher order Finite Element Methods are not monotone). But for 
fully nonlinear or degenerate elliptic, the only convergence proof currently available 
requires monotone schemes. One way to ensure monotone schemes is to use Wide 
Stencil Finite difference schemes, this has been done for the equation for motion by 
mean curvature, [Obe04], for the Infinity Laplace equation [Obe05], for functions 
of the eigenvalues [Obe08b], for Hamilton- Jacobi-Bellman equations [BZ03], and 
for the convex envelope [Obe08a]. Even for linear elliptic equations, a wide stencil 
scheme maybe necessary for to build a monotone scheme [MW53]. In some cases, 
simple finite difference schemes, with minor medications, can give good results, as is 
the case for the Monge- Ampere equation [BFOIO]. But we show below that simple 
finite difference schemes are not convergent for the Infinity Laplace equations. 

The second challenge, which is quite distinct from the first, is to build solvers 
for the finite difference schemes. For fixed values of /i, the fully the finite difference 
scheme is a finite dimensional nonlinear algebraic equation which must be solved. 
Build solvers demands very different techniques, and little progress has been made, 
in part, due to the fact that the discrete equations can be non- different iable, which 
precludes the use of the Newton's method. To date, the only general solver available 
is a fixed point iteration, which corresponds to solving the parabolic version of the 
equation for long time. This method is restricted by a nonlinear version of the 
CFL condition [Obe06] , which means the number of iterations required to solve the 
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equation increases with the problem size. For the Monge- Ampere equation, fast 
solvers have been built using Newton's method [FOllb] [FOlla], but this equation 
has a different structure (convex, different iable) from the Infinity Laplace (or the 
p-Laplace) equation. 

1.2. Contribution of this work. The first contribution of this work is to build a 
provably convergent discretization of the operator. The issue here is to ensure that 
the discretization convergences (in the limit of the discretization parameters going 
to zero) to the unique viscosity solution of the PDF. Simply using standard finite 
differences fails to converge, as shown below. 

The appropriate notion of weak solutions for the PDF is provided by viscosity 
solutions [CIL92, ACJ04]. The only schemes which can be proven to converge to vis- 
cosity solutions are monotone schemes [BS91]; these schemes satisfy the maximum 
principle at the discrete level [Obe06]. For the variational p-Laplacian, Galerkin 
Finite Flement methods could be used. But the game-theoretical version is not 
a divergence structure operator, so there is not a natural version of weak so- 
lutions. Monotone schemes can be proven to converge for the game-theoretical 
p-Laplacian (pLap). We prove convergence of the solution of the Wide Stencil 
finite difference schemes to the unique viscosity solution of the underlying equa- 
tion (pLap) (D). 

The second contribution of this work is to build fast solvers for (pLap). There 
are two reasonable ways to quantify the notion of a fast solver. The first notion 
of speed is absolute: the number of operations to solve the equation the should 
be proportional to the problem size. The second notion of speed is relative: we 
compare the speed of our solvers to the speed of solvers for a related but easier 
problem. Here, it is natural to compare with the solution speed of the Laplace 
equation. 

Fxplicit solvers are available and simple to implement, but they are not fast. Any 
monotone scheme can be solved using an iterative, explicit method [Obe06]. The 
explicit method can be interpreted as a Gauss-Seidel solver, or the forward Fuler 
method for the equation Ut = Ap. However the time step for the Fuler method is 
0{h?)^ where h is the spatial resolution. The explicit method is not fast because 
the number of iterations required for it to converge is 0{l/h^)^ which increases 
with the problem size. 

The method we propose is semi-implicit, with the implicit step given by solving 
the Laplace equation. The Laplace equation can be solved in 0{N) operations, 
using Fast Fourier Transforms, or 0{N\ogN)^ using sparse linear algebra. Thus, 
to the extent needed for our rather coarse analysis, both notions of speed coincide, 
provided the solution is obtained in a finite (small) number of iterations. 

1.3. The setting for the PDE. This work is concerned with the efficient numer- 
ical solution of a nonlinear, degenerate elliptic Partial Differential Fquation (PDF), 
the normalized Infinity Laplacian. The PDF operator is given by 
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(IL) Aoo^ = p^^;^ ^ U^,U:,^x,U^ 

where u{x) : R'^ ^ R. 
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We also study a closely related PDE, the game theoretical p-Laplacian, which 
interpolates between the 1-Laplacian, Ai, 

(MC) Aiu = \\/u\ dW{\/u/\\/u\) = Au- A^oU 

and the infinity Laplacian, Aoo- Expanding the Ai operator above leads to the 
identity 

(1) A = Aoo + Ai, 

which we record for future use. The p-Laplacian is the p weighted average of the 
1- and oo-Laplacians, 

(pLap) Ap = -Ai + -Aoo, P~^ + q~^ = 1- 

p q 

Special cases occur for p = 1 and oo, as above, and for p = 2 we obtain A2 = ^A. 
Using the identity (1) we can also write 

(2) Ap = aA^ /3Aoo, a = 1/p, (3 = {p - 2)/p. 

We consider the Dirichlet problem for the operator, in a domain (^ C M'^, with a 
given right hand side function g. 

(PDE) Apu{x) = g{x), for x e n, 

(D) u{x) = h{x)^ for X G dft. 

Here g represents a running cost for a probabilistic game [PS08]. When ^ = 0, the 
operator coincides with the variational p-Laplacian. The relationship between the 
game theoretical p-Laplacian and the variational p-Laplacian is given below. 

1.4. Interpretation of the operators. Here we define the game theoretical p- 
Laplacian operators. We are following the definitions given in the for the proba- 
bilistic games interpretation [PS08] (see also [MPRIO]). The normalized versions 
of the operators are also used in image processing [CMS98], [CEPB04], [MST06]. 

The variational p-Laplacian, which differs from our definition of the p-Laplacian 
by a factor involving the norm of the gradient of u^ is given by, 

dlY {\\/u{x)\P-^\/u) . 

It arises in the variational problem 

J[u] = / \Vu{x)\^ dx^ u = h on dQ 
Jn 

for 1 < p < 00. The Euler- Lagrange equation for the minimizer is 

diY{\\/u{x)\P-^\/u) =0. 

When the right hand side function, ^, is zero, the solutions of the both versions of 
the p-Laplacian equations coincide. 
An equivalent definition to (pLap) is 

(3) A,u=^^^^div{\yuix)r'Vu) 

To check consistency with definition (pLap), expand the two terms in the operator 
above to give ApU = -Au + ^^~ A^qU which is the representation (2). 
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2. Viscosity solutions 



In this section we recall the definition of viscosity solutions for the infinity Lapla- 
cian and the definition of consistency used in the convergence theory. 

Definition 1 (Viscosity solutions). (1) A continuous function u defined on the set 
[/, is a viscosity suhsolution of — Aoo'^ = in /7, if for every local maximum point 
X G U of u — (j)^ where (j) is C^ in some neighborhood of x, we have 

f-Aoo^(x) <0, ifL)^(x)^0, 

1 —r]^D'^(j){x)r] < for some |7^| < 1, if D(j){x) = 0. 

(2) A continuous function u defined on the set /7, is a viscosity super solution of 
— Aoo^ = in [/, if for every local maximum point x G U of u — (j), where (j) is C^ 
in some neighborhood of x, we have 

f-Aoo(/>(:r) >0, ifD(l){x)^0, 

\-r]^D^(l){x)r] > for some \r]\ < 1, if D^(x) = 0. 

(3) Moreover, a continuous function defined on the set [/, is a viscosity solution of 
— Aoo'^ = in [/, if it is both a viscosity subsolution and a viscosity supersolution 
in U. 

Consistency requires only that we can apply the test function definition in the 
limit. 

Definition 2 (Consistency). The numerical scheme A^'^^ is consistent if for every 
(j) e C^iU), and for every x eU, 

lim Al-''^^(,^)(x) = -Aoo.^(x) 

dx,dO^O 

if D(f){x) ^ 0, and 

(4) A < liminf A^'^^((^)(x) < limsup A^'^^(0)(x) < A 

dx,de^o dx,de^o 

where A, A are the least and greatest eigenvalues of D'^(j){x)^ otherwise. 

By a theorem of Barles-Souganidis [BS91], consistent, monotone schemes con- 
verge to the unique viscosity solution of the PDE. 

3. Discretization 

In this section we present the discretization of the Infinity Laplace operator, 
which is needed for convergence to the viscosity solution. The discretization we 
present here is different from the one in [Obe05]. The previous scheme was given 
by solving the a discrete version of the Lipschitz extension problem. This scheme 
simpler, since it separates and it also gives the correct scaling in h which is needed 
for a non-zero right hand side. 

By now it is well known that, for smooth functions with non- vanishing gradient, 
the operator is approximated by the average below. This result follows from Taylor 
expansions, and the fact that the minimum (or maximum) is in the direction of 
the gradient, at least to 0{h). We show below that the accuracy is actually 0{h?)^ 
which is an improvement over previous results. 
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Lemma 1. Let u{x) be a smooth function with non- vanishing gradient at x. Then 

.. A / N . u(y) — u(x) u(y) — u(x) ^, 9, 

(5) A^u{x)= min ^^\ ^ ^ + max ^ \ ^ ^ ^0{e^). 

\y-x\=e e \y-x\=e 6 

Proof. We prove the result in two dimensions, which is ah that is needed here. 
A longer proof is possible which in higher dimensions, which requires a Lagrange 
multiplier, A for the constraint, and an asymptotic expansion in A as well. 

It is enough to consider u{x) = p^x + ^x^Qx^ for nonzero p = Vix(x), and 
symmetric matrix Q = D'^u{x). Use the notation p = &, P^ = (— P2,Pi), for 
P = (^17^2)- With this notation, the operator, (IL) is given by 

(6) /^oou = {pfQp 

Write 

^2 

F{0) = u{ex{0)) = ep'^x{0) + -x{OfQx{0), x{0) = (cos((9), sin(6>)) 

Then a critical point of F is given by 

= F\0) = (p + eQx{0)yx\0) = (p + eQx^x^, 

where x'{6) = x^{6) = (— sin(^), cos(^)). Perform an asymptotic expansion the 
condition in e, with x = xq -\- exi, to obtain 

{p + eQxo^ixo + exi)^ + (9(e^) = 0. 

Then the 0(1) terms give 



which yields 

and the 0{e) terms give 

which yields 



/x^ = 



Xo = ±p, 
p'^xi + {x^Q)'^xo = 

Xi = p^ 

\p\ 



(Note that we are violating the constraint that x be a unit vector, but the constraint 
is stih satisfied to 0{e^)). Write 

x~^ = argmaxF(^), x~ = argminF(^) 



Then we have 



X = ±p + ecp , c = - 



\p\ • 

Inserting the values for x^, into the expression on the right hand side of (5) gives 

u(y) — u(x) u(y) — u(x) F(ex~) -\- F(ex~^) 

min ^^^ , ^ ^ + max - ^^^ v ; _ v ; v j 



\y-x\=e e^ \y-x\=e e 



. , 2f ip^VQp \\ 

■AooU + e^i —. Aiw 



which gives the desired result. D 
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Figure 1. Stencils for the 5, 9, and 17 point schemes 



However, for any given grid, it is impossible to sample the values on the entire 
circle. Instead, only values in a discrete set of directions are available. While it 
is certainly possible to interpolate the values onto the ball, quadratic interpolation 
is not monotone, so it violates the maximum principle, which is needed for the 
convergence proof. As we show in an example below, non-monotone schemes do 
not converge for this equation. 

So an additional discretization parameter is needed, which we present in what 
follows. 

Definition 3 (Spatial and directional resolution). Given a stencil of neighbouring 
grid points vi^. . . ^Vn on a Cartesian grid, define the local spatial resolution, h, to 
be the maximum length of the neighbours 



(7) 



h 



max \Vi\, 



and the local directional resolution, dO, to be the maximum directional distance to 
a neighbour 



(8) 



de 



n 




Vi 


max mm 


V - 




\v\ = l i=l 




\^i\ 



The direction vectors used will be on a grid, arranged as in Figure 1. In practice, 
we obtain acceptable accuracy using a relatively narrow stencil. 

Definition 4 (Scheme definition). Define the discretization of Aoo to be given by 



(9) 



A^^'u{x) = max 



u{x + Vi) — u{x) 



u{x + Vi) — u{x) 



where {vi} are the neighbours of the point x, as in Figure 1. 

Next we prove consistency, when the parameters dx,dO go to 0. For the grid 
points, we will assume 



(10) 

(11) 



whenever Vi is a neighbouring grid point, —Vi is as well, 

max,; I V,; I 

= o(l) 



mm^ \Vi^ 
Theorem 2. Let u he a smooth function in a neighbourhood of x, then 



(12) 



.h49 



u{x) = Aoo^(x) + O{d0 + h) 
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Proof. It is enough to consider 

u{x) = p^x + -x^Qx 

, for p = Vix(x), and symmetric matrix Q = D'^u{x). Use the notation p = ^. 
First consider the case Du ^ 0. Define 

u{x + Vi) — u{x) 

. u(x + Vi) — u(x) 
V = arg mm . — 

i \Vi\'^ 

Using a Taylor expansion for u^ we have 

u{x + Vi) - u{x) _ ViP + C\Vi\'^ 

RP " RP 

so, as /i ^ 0, the maximum occurs at v~^ = argmax^'O^p which is the closest 
direction vector to p, so 

(13) ^=p^O{dO) 
by (8) and similarly 

(14) P =p^O{dO), 

In addition, since according to (10) we assume that the grid points are arranged 
symmetrically, we have 

(15) ^ = -P 

as /i, dO -^ 0. 

Then using the Taylor expansion for u, insert (13,14) into (9) to obtain 

^ h rift , ^ f v^ v~ \ Iv^Qv^ v~Qv~ ,^,,, 

^ ^ ^ Vl^ I l^"l / 2 |v+p l^-p ^ ^ 

lv+Qv+ l^-Qv- 
= 2^;^ + 2T^^^ ^ by (15) 

= pQp + 0{h + di9) by (13,14) 

which is consistent, and accurate to 0{h -\- dO). 

In the case Du{x) = 0, we are only required only to verify (4). In this case we 
can assume u{x) = \x^ Qx^ and compute 

u{x ^ Vi) — u{x) 1 ^T^^ 1. /,^/ 7^ 7 9n 

max ^ , '\^ ^^ = - maxi?/ QiTi = -A + OidO + h^) 

i pil^ 2 i 2 

and similarly 

. U(X -\- Vi) — U(X) 1 . ^T^^ l^ /^/ .^ . 9n 

min ^ -^ ^ = - minVi^Qvi = -A + (9(dl9 + /i^) 

i \Vi\^ 2 i 2 

where A, A are the smallest and largest eigenvalues of Q, respectively. So then 

A^J\{x) = ^{X + A) + 0{de + h^) 
which is consistent, according to (4). D 
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3.1. The discretization of the p-Laplacian. We wish to build a consistent, 
monotone scheme for 

Starting with the scheme for A^o, we can simply combine this with the standard 
finite differences for the Laplacian, 

(16) AV) = 4^9^, 

u{x, y) = J {u{x -^h,y) ^ u{x - h,y) ^ u{x, y ^ h) ^ u{x, y - h)) 

Then, using the characterization of monotone schemes from [Obe06] , the com- 
bined scheme is still monotone. 

Theorem 3 (Convergence). The solution of the difference scheme for the p- Laplacian, {2), 
which is given by a convex combination of the schemes for (IL)^ (9) and the standard 
finite difference scheme for the Laplacian, (16) 

converges (uniformly on compact sets) as h^dO ^ to the solution of (3). 

Proof The uniform convergence of solutions of consistent, monotone schemes fol- 
lows from the main result of [BS91], so by appealing to this result we need only 
establish consistency and montonicity of the schemes. 

Consistency of the discretization of the A^^^ operator follows from Theorem 2. 
Consistency of the discretization for (3) follows since it is a convex combination of 
Aoo and the consistent discretization of the Laplacian (16). 

The definition (9) expresses the scheme A^^^ as a nondecreasing combination 
of differences between u{xi) — u{x)^ where Xi are neighbouring grids points to the 
reference point x. This characterizes elliptic schemes, which were proven to be 
monotone in [Obe06]. Monotonicity of the discretization of the Ap operator is a 
consequence of the fact that it is a convex combination of two elliptic schemes is 
also elliptic, which was also shown in [Obe06]. 

Together these results prove convergence. D 

4. Solvers 

When we discretize (3) we obtain a system of nonlinear equations. These equa- 
tions inherit contraction properties from (3), namely stability in the maximum 
norm. Our first goal is solve the equations, and our second goal is to solve them 
quickly, ideally in 0{N) iterations, where N is the number of variables in the dis- 
crete system of equations. 

The main issues to consider are stability of the iteration, and the convergence 
rate. For implicit or semi-implicit schemes, we also need to consider solvability of 
the operator. This generally requires that the implicit operator be linear. 

4.1. Explicit methods. There are no general methods available for solving non- 
divergence form nonlinear elliptic equations. For monotone schemes, explicit meth- 
ods can be used [Obe06]. The explicit method in the case of (3) is 

ii^+^ = ii^ + piApu"" - g) 

where, to ensure stability, the artificial time step p is restricted to be the inverse 
of the Lipschitz constant of the scheme, regarded as a grid function. It can be 
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computed easily from the scheme; p — 0(h?). ExpUcit methods are slow because of 
this restriction, which can be regarded as a nonlinear CFL condition. The number 
of number of iterations required for convergence grows with the problem size. 
Fully implicit methods require the solution of nonlinear equations. In the case 

of (3), the fully implicit scheme =^^ = Ap['u"^+^] — ^ is unconditionally stable, 

but requires solving the equation 

^/^+^ - pApii^+^ = ix^ - pg. 

which has the same difficulties as solving the time-independent equation. 

4.2. Semi-implicit methods. We are motivated by the work [Sme03] which built 
a semi-implicit solver for the one-Laplacian, Ai, by treating the linear part of the 
operator implicitly. 

Write 

Aoo = -(A + (Aoo-Ai)) 
which follows from (1). This leads to the semi-implicit scheme 

(17) - Am"+i = -(Ai - Aoo)m" - 2g. 

The general case for Ap follows. 

Remark 1. In practice, we only use the discretization of Aoo and the identity (1), 
so no discretization of Ai is needed. 

4.3. The semi- implicit scheme for Ap. Rewrite the equation (pLap) symmet- 
rically in terms of p, q as follows 

Ap = ^"'^'^"' (Ai + Aoo) + " ^^ (Ai - Aoo), 
which gives 

Ap = -(A + /3(Ai-Aoo; 
The resulting scheme, which generalizes (17), is given by 

(Iteration) - A^i^+i = /3(Aoo - Ai)ii^ -2g, f3 = l--. 

P 

The convergence of the iteration depends on whether the operator 

/?(-A)-i(Aoo-Ai) 

is a contraction is some norm. Clearly we can focus on the case |/3| = 1 since the 
case 1/3 1 < 1 is easier. 

4.4. Proof of contraction for a linear model. It would be desirable to prove 
that the iteration defined by (Iteration) converges to the solution. But to do so 
requires proving it is a contraction in some norm. The operator involved is not 
monotone, so we are unable to prove it is a contraction in the uniform norm. 
The nonlinearity makes proving convergence in other norms difficult. However we 
present a heuristic for why the operator A~^(Aoo — Ai) is a contraction based on 
an analogy with a linear operator. Instead of the operator above, we'll perform the 
analysis for 

V ^xx ^yy ) \ ^xx ~r '^yy ) 
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which is a hnear operator which is also a difference of degenerate ehiptic operators. 
To be concrete, consider the case of the unit square in two dimensions. Use a 
simultaneous eigenfunction expansion, 

i 

M(t)i = \^^(j)i, where M = d^x 
N(j)i = y^(j)i^ where N = dyy 



with 



Then 



So 



M-N M - TV ^ ^ x2 j,2 



M^N M^N^^ '^' ^^ 'A? 

In fact, the contraction rate is given by 






max — -^ ^ < 1, on a finite grid 

■ ■ \2 \_ -,,2 ' <^ 



J 

So in this simplified linear setting, the operator is a contraction. 

5. Numerical Results 

5.1. Failure of the standard finite difference scheme. Here we motivate the 
need for a convergent scheme, by showing that the standard finite difference scheme 
fails to converge. 

A natural scheme is given by standard finite differences, along with a small 
regularization for the norm of the gradient. For this we use standard centered 
finite differences for Uxx^Uyy^Ux^Uy^ and the with the symmetric scheme for Uxy, 

Ux{x, y) =— - {u{x ^h,y) - u{x - h, y)) + 0{h^), 

Uxx{x, y) =7^ {u{x ^h,y) - 2u{x, y) + u{x - h, y)) + 0{h^), 

1 
Ah? 

- 77^ (^(^ -h,y^h)^u{x^h,y- h)) + 0{h^) 

and similarly for the Uy^ Uyy terms. In order to regularize the gradient, we replaced 
||Vix|p with max{/iM|V^i|p}. 

We computed the solution with boundary conditions corresponding to the exact 
Aronsson solution [Aro68]. The finite difference scheme presented above fails to 
converge, see Figure 2. In this case, the solution has the form \x\ — \y\ in the centre. 
In fact, it can be shown using symmetry considerations that \x\ — \y\ is an exact 
solution of the symmetric finite difference scheme. On the fiat parts, the operator 
is zero, so we only need to check the corners. In fact, u{x^ y) = \x\ and u{x^ y) = \y\ 
are also exact solutions. While other discretizations are possible which break this 
symmetry, we tried several other simple consistent finite difference schemes and 
were always able to find examples where they failed to converge. 



ixy{x, 2/) = + 777 {u{x -^h,y^h)^u{x-h,y- h)) 
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Figure 2. 



Standard finite difference solutions of the oo- 



Laplacian. Boundary data is a |x 



4/3 



\y 



^/^, the computed so- 



lution is incorrect, with a singularity of the form \x\ 
origin, (a) Surface Plot, (b) plot of 2i(x,0). 



\y\ at the 



5.2. Plots of solutions of Aoo with varying boundary data. In Figure 3 we 
plot solutions of the oo-Laplacian, with variable boundary data, given by 

, , , , . 3x + 2?/ 



1^1 



8V14 



for 



0,1,2,3. 



Displayed are surface plots of the solution and a contour plot of the norm of the 
gradient of the solutions. Note how the changing boundary data moves the kinks in 
the solution from a symmetric arrangement to a perturbed one. The computational 
domain was given by [—1, 1]^ with a 200^ grid. 

5.3. Plots of solutions for varying p. We solved the p-Laplacian on a 200^ grid, 
with values of /3 = 0, 1/3, 2/3.1, which corresponds to p = 2, 3, 6, oo. The boundary 
data was \x\ — \y\ on [—1, 1]^ . 

In Figure 4 we show the surface plots of the solution and a contour plot of the 
norm of the gradient of the solution. The contours plotted are at levels sets equally 
spaced between and 1. Note in the surface plot, the sharpening of the saddle to 
a kink in the gradient as p increases. In the contour plots, the contours get smaller 
as p increases, showing that the solution is decreasing the size of the gradient. Note 
also that the shape of the contours go from circular (for the Laplacian) to diamond 
shaped (for the Infinity Laplacian). Also note that the gradient in larger along the 
axes, which is where (presumably) the singularity in the solution is found. 



6. Solver speed 

In this section we report on computational experiments which demonstrate the 
performance of the semi-implicit method. 

The first section compares the semi-implicit method, which convergences at a 
rate which is independent of the problems, to the explicit method, which takes 
increasingly longer to converge as the problem size grow. 

The next section gives further details on the solution of Aoo- Engineering preci- 
sion is obtained in the first few iterations. But numerical precision can still require 
many iterations. 
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Figure 3. Solutions of the oo-Laplacian. Boundary data is a 
\x\ - \y\ + Ci(3a; + 2y)/\/l4, for a = 0, 1/8, 2/8, 3/8, 4/8. 



FINITE DIFFERENCE METHODS FOR IFINITY LAPLACE EQUATION 



15 









5 0.5 -0.5 0.5 





5 0.5 -0.5 0.5 



Figure 4. Solution of the p-Laplacian with b = 0,1/3,2/3.1, 
giving p = 2,3,6, oo. Boundary data is \x\ — \y\. (a) Surface plot 
of the solutions, (b) Contour plots of the norm of gradient of the 
solution. 
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The final section compares the speed of the solution of the p-Laplacian. Setting 
p < oo improves the convergence rate of the solver. 

6.1. Speed of the explicit method. The explicit method converges exponen- 
tially, but at a rate which becomes depends on the problem size. For small problem 
sizes, n = 64, we get an error of .1 after 200 iterations, and .001 after 1000 it- 
erations. But for large problems sizes, n — 512, the error is .7 even after 1000 
iterations. This is prohibitively slow for larger problems sizes. 

In Figure 5 the error is shown as a function of the number of iterations, for 
different problem sizes. Increasing the problem size means the error increases as 
well, for a fixed number of iterations. 

We compared the semi-implicit method with the explicit method. The error 
decreases faster for the semi-implicit method. See Figure 6 for a typical example. 




Figure 5. Maximum error as a function of the number of it- 
erations in the solution of Aoo using the explicit method. For 
n = 64, 128, 256, 512. Errors are increasing with n. 




Figure 6. Maximum error as a function of iterations for solutions 
of Infinity Laplace. Semi-implicit method (better accuracy) and 



Explicit method (less accurate), with n 

,4/3 



128. Solution is x'^l'^ 



y 
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6.2. Solver speed for Infinity Laplace. In the section we examine the perfor- 
mance of the semi-imphcit method for solutions of Infinity Laplace. 

The advantage of the semi-implicit solver is that the convergence rate does not 
depend on the size of the problem. Another advantage is that the semi-implicit 
solver reduces the error to engineering tolerances, .1, after one iteration, and 
achieves close to another order of magnitude after four iterations. See Figure 7. 
This is useful for applications, for example in imaging, where high accuracy is not 
required. 




Figure 7. Maximum error as a function of the number of itera- 
tions in the solution of A^o using the semi-implicit method. For 
n = 64, 128, 256, 512. Errors are nearly independent of n. 



However, the specific rate of convergence varies for different solutions. Compar- 



ing several solutions, it appears to be slowest for the solution u{x^ y) 
see Figure 8 



,4/3 _ y4/3 



6.3. Solver speed for solutions of p-Laplacian. We study the convergence 
rate of the iterative method for the p-Laplacian with p < oo. We also study the 
convergence rate in terms of problem size. One notable fact is that the accuracy 
after just a few iterations is to engineering tolerances. 

First we display the speed of convergence in terms of the number of iterations 
for different values of a, where a = 1/p as given in (2). 

The results which study the convergence rate for a particular example solution as 
a function of a is presented in Figure 9 and Figure 10. The convergence rate is in- 
dependent of the size of the problem. The convergence increases with a, which is to 
be expected, since the operator Ap is becoming better approximated by the Lapla- 
cian. The convergence rate is consistent with an exponential error proportional to 
IQ^j'MN ^ where the convergence rate /i(a) is consistent with an affine expression 

/i(a) = -(ci +C2a) 

where ci is small and and C2 ~ 1. 

7. Conclusions 

We built convergent discretizations and fast solvers for the Aoo and A^ operators. 

Since solutions of the PDF (3) can be singular, it is important to use a convergent 

discretization. The theory of viscosity solutions provides the appropriate notion of 
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Figure 8. Maximum error as a function of iterations for solutions 
of Infinity Laplace. For n = 64, 128, 256. Semi-implicit method, 
using exact solution of (a) x^ — y'^ . (b) x^/^ — ?/^/^ (c) \x\ — \y\ 



Convergence 




40 60 

Iterations 



Figure 9. Convergence in the maximum norm as a function of 
the number of iterations, for a = 1/2, 1/2^, . . . , 1/2^^, 0. 



weak solutions for the PDF (3). We simplified and generalized the Wide Stencil 
finite difference scheme for Aoo used in [Obe05] to build a monotone finite difference 
scheme for Ap. The discretization took advantage of the fact the the Ap operator 
is a positive combination of A and Aoo, given by (2). 

A convergent discretization provides a finite dimensional equation whose so- 
lutions approximate the solutions of the PDF. Proving convergence is established 
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0.005 0.01 0.015 0.02 0.025 0.03 0.035 



Figure 10. Convergence rate as a function of a for a 
1/2, 1/2^, . . . , l/2l^ 0. The line of best fit is -.87a - .0085. 



using PDE techniques, which does not address efficiency of solution methods. Stan- 
dard explicit methods are slow, in the sense that the number of iterations required 
to approximately solve the equations to a given level of precision grows with the 
problem size. We introduced a semi-implicit solver, which converged at a rate inde- 
pendent of the problem size. The idea for the solver is to approximate the operator 
Ap by the Laplacian, and solve a linear equation at each iteration. While each it- 
eration is more costly than the explicit method, the added cost is more than made 
up for by the fact that accurate solutions can be obtained in a few iterations. 

Numerical results show that the method converges exponentially, at a rate which 
is independent of the problem size, even in the case p = oo. However, the conver- 
gence rate increases approximately linearly with 1/p, which is to be expected, since 
as p ^ 2 we recover the Laplacian and the solution is obtained in one step. For 
Of = or nearby, engineering precision (maximum error of 10~^) is obtained in a 
few iterations, but the convergence to numerical precision is slower. For a away 
from 0, the method is very fast, with numerical precision in about 20 iterations. 
The convergence rate also depends on the particular solution. 
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